Analysis for SSFA project: Two factor model

Filtered strongly by < 5% NMP

Table of contents

1. Used packages
2. Global settings
3. Load data
4. Model specification
5. Inference
   1. epLsar
   2. R²
   3. Asfc
   4. Smfc
   5. HAsfc9
   6. HAsfc81
6. Summary

Used packages

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
import pickle
import arviz as az
import pymc3 as pm

from matplotlib.colors import to_rgb

import scipy.stats as stats 

from IPython.display import display

import matplotlib as mpl

%load_ext autoreload
%autoreload 2

import plotting_lib

Global settings

Output

writeOut = True
outPathPlots = "../plots/statistical_model_two_factors_filter_strong/"
outPathData = "../derived_data/statistical_model_two_factors_filter_strong/"

Plotting

widthMM = 190 
widthInch = widthMM / 25.4
ratio = 0.66666
heigthInch = ratio*widthInch

SMALL_SIZE = 8
MEDIUM_SIZE = 10
BIGGER_SIZE = 12

plt.rc('font', size=SMALL_SIZE)          # controls default text sizes
plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title
plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels
plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize
plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title
sns.set_style("ticks")

dpi = 300

sizes = [SMALL_SIZE,MEDIUM_SIZE,BIGGER_SIZE]

Computing

numSamples = 1000
numCores = 4
numTune = 1000
numPredSamples = 2000
random_seed=36534535
target_accept = 0.99

Load data

datafile = "../derived_data/preprocessing/preprocessed_filter_strong.dat"

with open(datafile, "rb") as f:
    x1,x2,_,df,dataZ,dictMeanStd,dictTreatment,dictSoftware = pickle.load(f)    

Show that everything is correct:

display(pd.DataFrame.from_dict({'x1':x1,'x2':x2}))

	x1	x2
0	0	5
1	1	5
2	0	5
3	1	5
4	0	5
...	...	...
225	1	9
226	0	9
227	1	9
228	0	9
229	1	9

230 rows × 2 columns

x1 indicates the software used, x2 indicates the treatment applied.

for surfaceParam,(mean,std) in dictMeanStd.items():
    print("Surface parameter {} has mean {} and standard deviation {}".format(surfaceParam,mean,std))

Surface parameter epLsar has mean 0.0033908522378621737 and standard deviation 0.001989767613698366
Surface parameter Rsquared has mean 0.9973621855273914 and standard deviation 0.007969957191822462
Surface parameter Asfc has mean 13.800529685456086 and standard deviation 12.970336908999228
Surface parameter Smfc has mean 1.350346875497435 and standard deviation 7.832927757231397
Surface parameter HAsfc9 has mean 0.4835258022408901 and standard deviation 0.8643524891543747
Surface parameter HAsfc81 has mean 1.0354777484507558 and standard deviation 2.586788670504154

for k,v in sorted(dictTreatment.items(), key=lambda x: x[0]):    
    print("Number {} encodes treatment {}".format(k,v))

Number 0 encodes treatment BrushDirt
Number 1 encodes treatment BrushNoDirt
Number 2 encodes treatment Clover
Number 3 encodes treatment Clover+dust
Number 4 encodes treatment Control
Number 5 encodes treatment Dry bamboo
Number 6 encodes treatment Dry grass
Number 7 encodes treatment Dry lucerne
Number 8 encodes treatment Grass
Number 9 encodes treatment Grass+dust
Number 10 encodes treatment RubDirt

for k,v in dictSoftware.items():
    print("Number {} encodes software {}".format(k,v))

Number 0 encodes software ConfoMap
Number 1 encodes software Toothfrax

display(dataZ)

	index	TreatmentNumber	SoftwareNumber	DatasetNumber	NameNumber	epLsar_z	Rsquared_z	Asfc_z	Smfc_z	HAsfc9_z	HAsfc81_z
0	0	5	0	0	98	0.515753	0.080435	-0.165334	-0.134789	-0.402292	-0.259843
1	1	5	1	0	98	0.668494	0.275762	-0.268289	-0.157173	-0.402677	-0.260418
2	2	5	0	0	99	1.243847	-0.146130	0.128490	-0.134789	-0.360528	-0.285754
3	3	5	1	0	99	1.238912	0.264094	-0.046359	-0.157173	-0.363618	-0.299313
4	4	5	0	0	100	0.829617	-0.323025	-0.046215	-0.134789	-0.257176	-0.243567
...	...	...	...	...	...	...	...	...	...	...	...
225	225	9	1	2	52	0.519230	0.250668	-0.923702	0.007447	0.408438	0.029468
226	226	9	0	2	53	0.005954	0.186666	-0.843131	-0.110277	0.178138	0.194629
227	227	9	1	2	53	-0.126574	0.298222	-0.851345	-0.155367	0.546274	0.307811
228	228	9	0	2	54	-1.091196	0.265539	-0.949772	-0.014439	-0.130543	-0.089431
229	229	9	1	2	54	-1.350837	0.313027	-0.953590	-0.120251	-0.158980	-0.075589

230 rows × 11 columns

display(df)

	Dataset	Name	Software	Diet	Treatment	Before.after	NMP	NMP_cat	epLsar	Rsquared	Asfc	Smfc	HAsfc9	HAsfc81	NewEplsar	TreatmentNumber	SoftwareNumber	DatasetNumber	NameNumber
0	GuineaPigs	capor_2CC6B1_txP4_#1_1_100xL_1	ConfoMap	Dry bamboo	Dry bamboo	NaN	0.717312	0-5%	0.004417	0.998003	11.656095	0.294557	0.135803	0.363319	0.019460	5	0	0	98
1	GuineaPigs	capor_2CC6B1_txP4_#1_1_100xL_1	Toothfrax	Dry bamboo	Dry bamboo	NaN	0.717312	0-5%	0.004721	0.999560	10.320730	0.119219	0.135471	0.361833	NaN	5	1	0	98
2	GuineaPigs	capor_2CC6B1_txP4_#1_1_100xL_2	ConfoMap	Dry bamboo	Dry bamboo	NaN	1.674215	0-5%	0.005866	0.996198	15.467083	0.294557	0.171903	0.296292	0.020079	5	0	0	99
3	GuineaPigs	capor_2CC6B1_txP4_#1_1_100xL_2	Toothfrax	Dry bamboo	Dry bamboo	NaN	1.674215	0-5%	0.005856	0.999467	13.199232	0.119219	0.169232	0.261217	NaN	5	1	0	99
4	GuineaPigs	capor_2CC6B1_txP4_#1_1_100xL_3	ConfoMap	Dry bamboo	Dry bamboo	NaN	1.760409	0-5%	0.005042	0.994788	13.201101	0.294557	0.261235	0.405422	0.019722	5	0	0	100
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
225	Sheeps	L8-Ovis-90730-lm2sin-a	Toothfrax	Grass+dust	Grass+dust	NaN	0.000000	0-5%	0.004424	0.999360	1.819802	1.408678	0.836560	1.111706	NaN	9	1	2	52
226	Sheeps	L8-Ovis-90764-lm2sin-a	ConfoMap	Grass+dust	Grass+dust	NaN	0.000000	0-5%	0.003403	0.998850	2.864831	0.486556	0.637499	1.538943	0.018978	9	0	2	53
227	Sheeps	L8-Ovis-90764-lm2sin-a	Toothfrax	Grass+dust	Grass+dust	NaN	0.000000	0-5%	0.003139	0.999739	2.758297	0.133366	0.955699	1.831721	NaN	9	1	2	53
228	Sheeps	L8-Ovis-90814-lm2sin-a	ConfoMap	Grass+dust	Grass+dust	NaN	0.000000	0-5%	0.001220	0.999479	1.481662	1.237247	0.370691	0.804138	0.017498	9	0	2	54
229	Sheeps	L8-Ovis-90814-lm2sin-a	Toothfrax	Grass+dust	Grass+dust	NaN	0.000000	0-5%	0.000703	0.999857	1.432148	0.408433	0.346111	0.839946	NaN	9	1	2	54

230 rows × 19 columns

Model specification

class TwoFactorModel(pm.Model):
    
    """
    Compute params of priors and hyperpriors.
    """
    def getParams(self,x1,x2,y):
        # get lengths
        Nx1Lvl = np.unique(x1).size
        Nx2Lvl = np.unique(x2).size
                        
        dims = (Nx1Lvl, Nx2Lvl)
        
        ### get standard deviations
        
        # convert to pandas dataframe to use their logic
        df = pd.DataFrame.from_dict({'x1':x1,'x2':x2,'y':y})
        
        s1 = df.groupby('x1').std()['y'].max()
        s2 = df.groupby('x2').std()['y'].max()
                
        stdSingle = (s1, s2)
        
        prefac = 0.05
        s12 = prefac * np.linalg.norm([s1,s2])
        
        stdMulti = (s12)
        
        return (dims, stdSingle, stdMulti)
    
    def printParams(self,x1,x2,y):
        dims, stdSingle, stdMulti = self.getParams(x1,x2,y)
        Nx1Lvl, Nx2Lvl = dims
        s1, s2 = stdSingle
        s12 = stdMulti
        
        print("The number of levels of the x variables are {}".format(dims))
        print("The standard deviations used for the beta priors are  {}".format(stdSingle))
        print("The standard deviations used for the M12 priors are {}".format(stdMulti))
    
    def __init__(self,name,x1,x2,y,model=None):
        
        # call super's init first, passing model and name
        super().__init__(name, model)
        
        # get parameter of hyperpriors
        dims, stdSingle, stdMulti = self.getParams(x1,x2,y)
        Nx1Lvl, Nx2Lvl = dims
        s1, s2 = stdSingle
        s12 = stdMulti
        
        ### hyperpriors ### 
        # observation hyperpriors
        lamY = 1/30.
        muGamma = 0.5
        sigmaGamma = 2.
        
        # prediction hyperpriors
        sigma0 = pm.HalfNormal('sigma0',sd=1)
        sigma1 = pm.HalfNormal('sigma1',sd=s1, shape=Nx1Lvl)       
        sigma2 = pm.HalfNormal('sigma2',sd=s2, shape=Nx2Lvl)
        beta2 = (np.sqrt(6)*sigma2)/(np.pi)
        
        mu_b0 = pm.Normal('mu_b0', mu=0., sd=1)
        mu_b1 = pm.Normal('mu_b1', mu=0., sd=1, shape=Nx1Lvl)        
        mu_b2 = pm.Normal('mu_b2', mu=0., sd=1, shape=Nx2Lvl)       
        
        sigma12 = pm.HalfNormal('sigma12',sd=s12)
                               
        ### priors ### 
        # observation priors        
        nuY = pm.Exponential('nuY',lam=lamY)
        sigmaY = pm.Gamma('sigmaY',mu=muGamma, sigma=sigmaGamma)
        
        # prediction priors
        b0_dist = pm.Normal('b0_dist', mu=0, sd=1)
        b0 = pm.Deterministic("b0", mu_b0 + b0_dist * sigma0)
       
        b1_dist = pm.Normal('b1_dist', mu=0, sd=1)
        b1 = pm.Deterministic("b1", mu_b1 + b1_dist * sigma1)
                        
        b2_beta = pm.HalfNormal('b2_beta', sd=beta2, shape=Nx2Lvl)
        b2_dist = pm.Gumbel('b2_dist', mu=0, beta=1)
        b2 = pm.Deterministic("b2", mu_b2 + b2_beta * b2_dist)
        
        mu_M12 = pm.Normal('mu_M12', mu=0., sd=1, shape=[Nx1Lvl, Nx2Lvl])
        M12_dist = pm.Normal('M12_dist', mu=0, sd=1)
        M12 = pm.Deterministic("M12", mu_M12 + M12_dist * sigma12)
        
        #### prediction ###         
        mu = pm.Deterministic('mu',b0 + b1[x1]+ b2[x2] +  M12[x1,x2] )
                                        
        ### observation ### 
        y = pm.StudentT('y',nu = nuY, mu=mu, sd=sigmaY, observed=y)

Inference

epLsar

with pm.Model() as model:
    epLsarModel = TwoFactorModel('epLsar',x1,x2,dataZ.epLsar_z.values)

Verify model settings

epLsarModel.printParams(x1,x2,dataZ.epLsar_z.values)

The number of levels of the x variables are (2, 11)
The standard deviations used for the beta priors are  (1.022692176080412, 1.4502540614172976)
The standard deviations used for the M12 priors are 0.08872902751740064

try:
    graph_epLsar = pm.model_to_graphviz(epLsarModel)    
except:
    graph_epLsar = "Could not make graph"
graph_epLsar

Check prior choice

with epLsarModel as model:
    prior_pred_epLsar = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)

plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_epLsar,dataZ.epLsar_z.values,'epLsar')

Prior choice is as intended: Broad over the data range.

Sampling

with epLsarModel as model:
    trace_epLsar = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.95,random_seed=random_seed)
    #fit_epLsar = pm.fit(random_seed=random_seed)
    #trace_epLsar = fit_epLsar.sample(draws=numSamples)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [epLsar_M12_dist, epLsar_mu_M12, epLsar_b2_dist, epLsar_b2_beta, epLsar_b1_dist, epLsar_b0_dist, epLsar_sigmaY, epLsar_nuY, epLsar_sigma12, epLsar_mu_b2, epLsar_mu_b1, epLsar_mu_b0, epLsar_sigma2, epLsar_sigma1, epLsar_sigma0]

100.00% [8000/8000 01:29<00:00 Sampling 4 chains, 1 divergences]

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 90 seconds.
There was 1 divergence after tuning. Increase `target_accept` or reparameterize.

with epLsarModel as model:
    if writeOut:
        with open(outPathData + 'model_{}.pkl'.format('epLsar'), 'wb') as buff:
            pickle.dump({'model':epLsarModel, 'trace': trace_epLsar}, buff)            

Save for later comparison

if writeOut:
    np.save('../derived_data/statistical_model_two_factors_filter_strong/epLsar_oldb1', trace_epLsar['epLsar_b1'])
    np.save('../derived_data/statistical_model_two_factors_filter_strong/epLsar_oldb2', trace_epLsar['epLsar_b2'])
    np.save('../derived_data/statistical_model_two_factors_filter_strong/epLsar_oldM12', trace_epLsar['epLsar_M12'])

Check sampling

with epLsarModel as model:
    dataTrace_epLsar = az.from_pymc3(trace=trace_epLsar)

pm.summary(dataTrace_epLsar,hdi_prob=0.95).round(2)

	mean	sd	hdi_2.5%	hdi_97.5%	mcse_mean	mcse_sd	ess_mean	ess_sd	ess_bulk	ess_tail	r_hat
epLsar_mu_b0	-0.06	0.80	-1.62	1.54	0.01	0.01	6568.0	1939.0	6611.0	2872.0	1.0
epLsar_mu_b1[0]	0.01	0.70	-1.47	1.30	0.01	0.01	5782.0	1770.0	5782.0	2684.0	1.0
epLsar_mu_b1[1]	-0.07	0.70	-1.52	1.25	0.01	0.01	6014.0	2052.0	5972.0	3415.0	1.0
epLsar_mu_b2[0]	-0.13	0.69	-1.52	1.18	0.01	0.01	5704.0	1652.0	5713.0	2643.0	1.0
epLsar_mu_b2[1]	-0.11	0.68	-1.42	1.25	0.01	0.01	6099.0	1286.0	6033.0	2908.0	1.0
...	...	...	...	...	...	...	...	...	...	...	...
epLsar_mu[225]	0.37	0.29	-0.20	0.95	0.00	0.00	4203.0	3833.0	4195.0	3730.0	1.0
epLsar_mu[226]	0.36	0.30	-0.25	0.94	0.00	0.00	4139.0	3594.0	4122.0	3705.0	1.0
epLsar_mu[227]	0.37	0.29	-0.20	0.95	0.00	0.00	4203.0	3833.0	4195.0	3730.0	1.0
epLsar_mu[228]	0.36	0.30	-0.25	0.94	0.00	0.00	4139.0	3594.0	4122.0	3705.0	1.0
epLsar_mu[229]	0.37	0.29	-0.20	0.95	0.00	0.00	4203.0	3833.0	4195.0	3730.0	1.0

334 rows × 11 columns

plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_epLsar,dataTrace_epLsar,'epLsar')

/home/bob/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/arviz/plots/backends/matplotlib/pairplot.py:216: UserWarning: rcParams['plot.max_subplots'] (40) is smaller than the number of resulting pair plots with these variables, generating only a 8x8 grid
  UserWarning,

with epLsarModel as model:
    plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_epLsar,'epLsar')

with epLsarModel as model:
    plotting_lib.pm.energyplot(trace_epLsar)

Posterior predictive distribution

with epLsarModel as model:
    posterior_pred_epLsar = pm.sample_posterior_predictive(trace_epLsar,samples=numPredSamples,random_seed=random_seed)

/home/bob/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/pymc3/sampling.py:1708: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample
  "samples parameter is smaller than nchains times ndraws, some draws "

100.00% [2000/2000 00:01<00:00]

plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_epLsar,posterior_pred_epLsar,dataZ.epLsar_z.values,'epLsar')

Level plots

plotting_lib.plotLevels(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_epLsar,'epLsar',x1,x2)

plotting_lib.plotLevelsStd(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_epLsar,'epLsar',x1,x2)

Posterior and contrasts

df_hdi_epLsar = plotting_lib.plotTreatmentPosterior(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_epLsar,'epLsar',x1,x2)

df_hdi_epLsar

	Treatment_i	Treatment_j	hdi_ConfoMap_2.5%	hdi_ConfoMap_97.5%	isSignificant_on_ConfoMap	hdi_Toothfrax_2.5%	hdi_Toothfrax_97.5%	isSignificant_on_Toothfrax
0	Dry grass	Dry bamboo	-0.004643	-0.002188	True	-0.004397	-0.001964	True
1	Dry lucerne	Dry bamboo	-0.003704	-0.001872	True	-0.003929	-0.002080	True
2	Dry lucerne	Dry grass	-0.000554	0.001779	False	-0.001026	0.001433	False
3	BrushNoDirt	BrushDirt	-0.001433	0.002212	False	-0.002135	0.001267	False
4	Control	BrushDirt	-0.001220	0.002299	False	-0.000928	0.002773	False
5	Control	BrushNoDirt	-0.001621	0.001920	False	-0.000468	0.002948	False
6	RubDirt	BrushDirt	-0.000880	0.003012	False	-0.002004	0.001707	False
7	RubDirt	BrushNoDirt	-0.001245	0.002679	False	-0.001468	0.001930	False
8	RubDirt	Control	-0.001353	0.002504	False	-0.002719	0.001017	False
9	Clover+dust	Clover	-0.001200	0.001581	False	-0.000904	0.001890	False
10	Grass	Clover	-0.000091	0.003047	False	0.000486	0.003751	True
11	Grass	Clover+dust	-0.000281	0.002971	False	0.000105	0.003410	True
12	Grass+dust	Clover	0.000179	0.003309	True	0.000571	0.003553	True
13	Grass+dust	Clover+dust	0.000003	0.003144	True	0.000194	0.003165	True
14	Grass+dust	Grass	-0.001516	0.001911	False	-0.001688	0.001791	False

plotting_lib.plotTreatmentPosteriorDiff(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_epLsar,'epLsar',x1,x2)

if writeOut:
    df_hdi_epLsar.to_csv(outPathData+ 'hdi_{}.csv'.format('epLsar'))

R²

with pm.Model() as model:
    RsquaredModel = TwoFactorModel('Rsquared',x1,x2,dataZ["Rsquared_z"].values)

Verify model settings

RsquaredModel.printParams(x1,x2,dataZ["Rsquared_z"].values)

The number of levels of the x variables are (2, 11)
The standard deviations used for the beta priors are  (1.3690772711592387, 3.503507127699258)
The standard deviations used for the M12 priors are 0.1880753490508813

pm.model_to_graphviz(RsquaredModel)

Check prior choice

with RsquaredModel as model:
    prior_pred_Rsquared = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)

plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Rsquared,dataZ["Rsquared_z"].values,'Rsquared')

Sampling

with RsquaredModel as model:
    trace_Rsquared = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.99,random_seed=random_seed)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [Rsquared_M12_dist, Rsquared_mu_M12, Rsquared_b2_dist, Rsquared_b2_beta, Rsquared_b1_dist, Rsquared_b0_dist, Rsquared_sigmaY, Rsquared_nuY, Rsquared_sigma12, Rsquared_mu_b2, Rsquared_mu_b1, Rsquared_mu_b0, Rsquared_sigma2, Rsquared_sigma1, Rsquared_sigma0]

100.00% [8000/8000 34:38<00:00 Sampling 4 chains, 0 divergences]

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 2079 seconds.

with RsquaredModel as model:
    if writeOut:
        with open(outPathData + 'model_{}.pkl'.format('Rsquared'), 'wb') as buff:
            pickle.dump({'model': RsquaredModel, 'trace': trace_Rsquared}, buff)

Check sampling

with RsquaredModel as model:
    dataTrace_Rsquared = az.from_pymc3(trace=trace_Rsquared)

pm.summary(dataTrace_Rsquared,hdi_prob=0.95).round(2)

	mean	sd	hdi_2.5%	hdi_97.5%	mcse_mean	mcse_sd	ess_mean	ess_sd	ess_bulk	ess_tail	r_hat
Rsquared_mu_b0	0.07	0.83	-1.56	1.67	0.01	0.01	4417.0	2226.0	4420.0	3217.0	1.0
Rsquared_mu_b1[0]	-0.04	0.73	-1.42	1.47	0.01	0.01	3480.0	2287.0	3469.0	2992.0	1.0
Rsquared_mu_b1[1]	0.09	0.72	-1.23	1.61	0.01	0.01	3364.0	2286.0	3364.0	2820.0	1.0
Rsquared_mu_b2[0]	0.03	0.66	-1.27	1.32	0.01	0.01	2953.0	2039.0	2997.0	2879.0	1.0
Rsquared_mu_b2[1]	0.03	0.67	-1.26	1.30	0.01	0.01	3269.0	2311.0	3301.0	2698.0	1.0
...	...	...	...	...	...	...	...	...	...	...	...
Rsquared_mu[225]	0.28	0.01	0.25	0.31	0.00	0.00	3553.0	3530.0	3554.0	3929.0	1.0
Rsquared_mu[226]	0.20	0.02	0.15	0.25	0.00	0.00	3916.0	3916.0	3926.0	3480.0	1.0
Rsquared_mu[227]	0.28	0.01	0.25	0.31	0.00	0.00	3553.0	3530.0	3554.0	3929.0	1.0
Rsquared_mu[228]	0.20	0.02	0.15	0.25	0.00	0.00	3916.0	3916.0	3926.0	3480.0	1.0
Rsquared_mu[229]	0.28	0.01	0.25	0.31	0.00	0.00	3553.0	3530.0	3554.0	3929.0	1.0

334 rows × 11 columns

plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_Rsquared,dataTrace_Rsquared,'Rsquared')

/home/bob/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/arviz/plots/backends/matplotlib/pairplot.py:216: UserWarning: rcParams['plot.max_subplots'] (40) is smaller than the number of resulting pair plots with these variables, generating only a 8x8 grid
  UserWarning,

with RsquaredModel as model:
    plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_Rsquared,'Rsquared')

with RsquaredModel as model:
    plotting_lib.pm.energyplot(trace_Rsquared)

Posterior predictive distribution

with RsquaredModel as model:
    posterior_pred_Rsquared = pm.sample_posterior_predictive(trace_Rsquared,samples=numPredSamples,random_seed=random_seed)

/home/bob/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/pymc3/sampling.py:1708: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample
  "samples parameter is smaller than nchains times ndraws, some draws "

100.00% [2000/2000 00:01<00:00]

plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Rsquared,posterior_pred_Rsquared,dataZ["Rsquared_z"].values,'Rsquared')

Compare prior and posterior for model parameters

with RsquaredModel as model:
    pm_data_Rsquared = az.from_pymc3(trace=trace_Rsquared,prior=prior_pred_Rsquared,posterior_predictive=posterior_pred_Rsquared)

arviz.data.io_pymc3 - WARNING - posterior predictive variable Rsquared_y's shape not compatible with number of chains and draws. This can mean that some draws or even whole chains are not represented.

plotting_lib.plotPriorPosteriorB(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,pm_data_Rsquared,'Rsquared')

plotting_lib.plotLevels(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_Rsquared,'Rsquared',x1,x2)

plotting_lib.plotLevelsStd(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_Rsquared,'Rsquared',x1,x2)

Posterior and contrasts

plotting_lib.plotPosterior(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,pm_data_Rsquared,'Rsquared')

df_hdi_R = plotting_lib.plotTreatmentPosterior(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_Rsquared,'Rsquared',x1,x2)

df_hdi_R

	Treatment_i	Treatment_j	hdi_ConfoMap_2.5%	hdi_ConfoMap_97.5%	isSignificant_on_ConfoMap	hdi_Toothfrax_2.5%	hdi_Toothfrax_97.5%	isSignificant_on_Toothfrax
0	Dry grass	Dry bamboo	0.001541	0.002830	True	-0.000339	0.000455	False
1	Dry lucerne	Dry bamboo	0.000551	0.001855	True	-0.000373	0.000178	False
2	Dry lucerne	Dry grass	-0.001538	-0.000439	True	-0.000609	0.000202	False
3	BrushNoDirt	BrushDirt	-0.001141	0.000105	False	-0.000936	0.000088	False
4	Control	BrushDirt	-0.001287	0.000661	False	-0.000664	0.000164	False
5	Control	BrushNoDirt	-0.000739	0.001189	False	-0.000262	0.000716	False
6	RubDirt	BrushDirt	-0.001555	-0.000071	True	-0.001123	0.000027	False
7	RubDirt	BrushNoDirt	-0.001066	0.000435	False	-0.000755	0.000524	False
8	RubDirt	Control	-0.001689	0.000501	False	-0.000874	0.000229	False
9	Clover+dust	Clover	-0.000440	0.000460	False	-0.000216	0.000461	False
10	Grass	Clover	-0.000685	0.000545	False	-0.000706	0.000045	False
11	Grass	Clover+dust	-0.000705	0.000535	False	-0.000831	-0.000136	True
12	Grass+dust	Clover	-0.000430	0.000526	False	-0.000451	0.000218	False
13	Grass+dust	Clover+dust	-0.000516	0.000495	False	-0.000568	0.000045	False
14	Grass+dust	Grass	-0.000502	0.000785	False	-0.000141	0.000578	False

plotting_lib.plotTreatmentPosteriorDiff(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_Rsquared,'Rsquared',x1,x2)

if writeOut:
    df_hdi_R.to_csv(outPathData+ 'hdi_{}.csv'.format('Rsquared'))

Asfc

with pm.Model() as model:
    AsfcModel = TwoFactorModel('Asfc',x1,x2,dataZ["Asfc_z"].values)

Verify model settings

AsfcModel.printParams(x1,x2,dataZ["Asfc_z"].values)

The number of levels of the x variables are (2, 11)
The standard deviations used for the beta priors are  (1.0817297140638054, 1.6678950572969222)
The standard deviations used for the M12 priors are 0.09939835381489054

pm.model_to_graphviz(AsfcModel)

Check prior choice

with AsfcModel as model:
    prior_pred_Asfc = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)

plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Asfc,dataZ["Asfc_z"].values,'Asfc')

Prior choice is as intended: Broad over the data range.

Sampling

with AsfcModel as model:
    trace_Asfc = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.99,random_seed=random_seed)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [Asfc_M12_dist, Asfc_mu_M12, Asfc_b2_dist, Asfc_b2_beta, Asfc_b1_dist, Asfc_b0_dist, Asfc_sigmaY, Asfc_nuY, Asfc_sigma12, Asfc_mu_b2, Asfc_mu_b1, Asfc_mu_b0, Asfc_sigma2, Asfc_sigma1, Asfc_sigma0]

100.00% [8000/8000 10:44<00:00 Sampling 4 chains, 0 divergences]

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 645 seconds.
The estimated number of effective samples is smaller than 200 for some parameters.

with AsfcModel as model:
    if writeOut:
        with open(outPathData + 'model_{}.pkl'.format('Asfc'), 'wb') as buff:
            pickle.dump({'model': AsfcModel, 'trace': trace_Asfc}, buff)

Check sampling

with AsfcModel as model:
    dataTrace_Asfc = az.from_pymc3(trace=trace_Asfc)

pm.summary(dataTrace_Asfc,hdi_prob=0.95).round(2)

	mean	sd	hdi_2.5%	hdi_97.5%	mcse_mean	mcse_sd	ess_mean	ess_sd	ess_bulk	ess_tail	r_hat
Asfc_mu_b0	-0.03	0.82	-1.62	1.50	0.01	0.01	4569.0	1742.0	4604.0	3062.0	1.0
Asfc_mu_b1[0]	0.06	0.72	-1.37	1.47	0.01	0.01	5112.0	1836.0	5115.0	2795.0	1.0
Asfc_mu_b1[1]	-0.11	0.72	-1.55	1.26	0.01	0.01	4876.0	2061.0	4878.0	2946.0	1.0
Asfc_mu_b2[0]	1.16	0.76	-0.36	2.66	0.01	0.01	3730.0	3282.0	3758.0	3092.0	1.0
Asfc_mu_b2[1]	0.90	0.75	-0.54	2.40	0.01	0.01	4520.0	3394.0	4559.0	2919.0	1.0
...	...	...	...	...	...	...	...	...	...	...	...
Asfc_mu[225]	-0.93	0.04	-1.02	-0.85	0.00	0.00	4032.0	4032.0	4049.0	3672.0	1.0
Asfc_mu[226]	-0.93	0.04	-1.01	-0.84	0.00	0.00	3668.0	3651.0	3647.0	3391.0	1.0
Asfc_mu[227]	-0.93	0.04	-1.02	-0.85	0.00	0.00	4032.0	4032.0	4049.0	3672.0	1.0
Asfc_mu[228]	-0.93	0.04	-1.01	-0.84	0.00	0.00	3668.0	3651.0	3647.0	3391.0	1.0
Asfc_mu[229]	-0.93	0.04	-1.02	-0.85	0.00	0.00	4032.0	4032.0	4049.0	3672.0	1.0

334 rows × 11 columns

plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_Asfc,dataTrace_Asfc,'Asfc')

/home/bob/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/arviz/plots/backends/matplotlib/pairplot.py:216: UserWarning: rcParams['plot.max_subplots'] (40) is smaller than the number of resulting pair plots with these variables, generating only a 8x8 grid
  UserWarning,

with AsfcModel as model:
    plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_Asfc,'Asfc')

with AsfcModel as model:
    plotting_lib.pm.energyplot(trace_Asfc)

Posterior predictive distribution

with AsfcModel as model:
    posterior_pred_Asfc = pm.sample_posterior_predictive(trace_Asfc,samples=numPredSamples,random_seed=random_seed)

/home/bob/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/pymc3/sampling.py:1708: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample
  "samples parameter is smaller than nchains times ndraws, some draws "

100.00% [2000/2000 00:01<00:00]

plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Asfc,posterior_pred_Asfc,dataZ["Asfc_z"].values,'Asfc')

plotting_lib.plotLevels(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_Asfc,'Asfc',x1,x2)

plotting_lib.plotLevelsStd(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_Asfc,'Asfc',x1,x2)

Compare prior and posterior for model parameters

with AsfcModel as model:
    pm_data_Asfc = az.from_pymc3(trace=trace_Asfc,prior=prior_pred_Asfc,posterior_predictive=posterior_pred_Asfc)

arviz.data.io_pymc3 - WARNING - posterior predictive variable Asfc_y's shape not compatible with number of chains and draws. This can mean that some draws or even whole chains are not represented.

plotting_lib.plotPriorPosteriorB(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,pm_data_Asfc,'Asfc')

Posterior and contrasts

plotting_lib.plotPosterior(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,pm_data_Asfc,'Asfc')

df_hdi_Asfc = plotting_lib.plotTreatmentPosterior(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_Asfc,'Asfc',x1,x2)

df_hdi_Asfc

	Treatment_i	Treatment_j	hdi_ConfoMap_2.5%	hdi_ConfoMap_97.5%	isSignificant_on_ConfoMap	hdi_Toothfrax_2.5%	hdi_Toothfrax_97.5%	isSignificant_on_Toothfrax
0	Dry grass	Dry bamboo	-2.914947	2.365610	False	-2.600433	1.769296	False
1	Dry lucerne	Dry bamboo	-1.399234	2.293610	False	-1.732294	1.659878	False
2	Dry lucerne	Dry grass	-1.688190	3.380153	False	-1.669210	2.232344	False
3	BrushNoDirt	BrushDirt	-15.196748	3.533049	False	-14.476359	3.598702	False
4	Control	BrushDirt	-35.751145	4.803498	False	-31.583528	3.896526	False
5	Control	BrushNoDirt	-33.267965	15.869328	False	-30.695311	14.827444	False
6	RubDirt	BrushDirt	-43.457907	-34.032025	True	-36.298904	-25.764600	True
7	RubDirt	BrushNoDirt	-43.669076	-23.369520	True	-36.904964	-15.465674	True
8	RubDirt	Control	-44.050859	-2.428273	True	-37.081175	0.973047	False
9	Clover+dust	Clover	-3.719954	2.410717	False	-3.596101	2.037584	False
10	Grass	Clover	-5.067046	1.295456	False	-4.849023	0.951618	False
11	Grass	Clover+dust	-3.219981	0.650326	False	-2.833752	0.633075	False
12	Grass+dust	Clover	-5.715432	0.022721	False	-5.235751	0.294740	False
13	Grass+dust	Clover+dust	-3.949583	-0.482832	True	-3.407486	-0.126538	True
14	Grass+dust	Grass	-2.676912	0.879903	False	-2.324307	1.051653	False

plotting_lib.plotTreatmentPosteriorDiff(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_Asfc,'Asfc',x1,x2)

if writeOut:
    df_hdi_Asfc.to_csv(outPathData+ 'hdi_{}.csv'.format('Asfc'))

Smfc

with pm.Model() as model:
    SmfcModel = TwoFactorModel('Smfc',x1,x2,dataZ.Smfc_z.values)

Verify model settings

SmfcModel.printParams(x1,x2,dataZ.Smfc_z.values)

The number of levels of the x variables are (2, 11)
The standard deviations used for the beta priors are  (1.2914215627674575, 3.238203806583398)
The standard deviations used for the M12 priors are 0.17431102622720077

pm.model_to_graphviz(SmfcModel)

Check prior choice

with SmfcModel as model:
    prior_pred_Smfc = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)

plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Smfc,dataZ.Smfc_z.values,'Smfc')

Prior choice is as intended: Broad over the data range.

Sampling

with SmfcModel as model:
    trace_Smfc = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.99,random_seed=random_seed)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [Smfc_M12_dist, Smfc_mu_M12, Smfc_b2_dist, Smfc_b2_beta, Smfc_b1_dist, Smfc_b0_dist, Smfc_sigmaY, Smfc_nuY, Smfc_sigma12, Smfc_mu_b2, Smfc_mu_b1, Smfc_mu_b0, Smfc_sigma2, Smfc_sigma1, Smfc_sigma0]

10.12% [810/8000 00:15<02:21 Sampling 4 chains, 0 divergences]

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
~/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/pymc3/sampling.py in _mp_sample(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, callback, discard_tuned_samples, mp_ctx, pickle_backend, **kwargs)
   1485             with sampler:
-> 1486                 for draw in sampler:
   1487                     trace = traces[draw.chain - chain]

~/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/pymc3/parallel_sampling.py in __iter__(self)
    491         while self._active:
--> 492             draw = ProcessAdapter.recv_draw(self._active)
    493             proc, is_last, draw, tuning, stats, warns = draw

~/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/pymc3/parallel_sampling.py in recv_draw(processes, timeout)
    351         pipes = [proc._msg_pipe for proc in processes]
--> 352         ready = multiprocessing.connection.wait(pipes)
    353         if not ready:

/usr/lib/python3.7/multiprocessing/connection.py in wait(object_list, timeout)
    919             while True:
--> 920                 ready = selector.select(timeout)
    921                 if ready:

/usr/lib/python3.7/selectors.py in select(self, timeout)
    414         try:
--> 415             fd_event_list = self._selector.poll(timeout)
    416         except InterruptedError:

KeyboardInterrupt: 

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
<ipython-input-121-4de91012caf3> in <module>
      1 with SmfcModel as model:
----> 2     trace_Smfc = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.99,random_seed=random_seed)

~/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/pymc3/sampling.py in sample(draws, step, init, n_init, start, trace, chain_idx, chains, cores, tune, progressbar, model, random_seed, discard_tuned_samples, compute_convergence_checks, callback, return_inferencedata, idata_kwargs, mp_ctx, pickle_backend, **kwargs)
    543         _print_step_hierarchy(step)
    544         try:
--> 545             trace = _mp_sample(**sample_args, **parallel_args)
    546         except pickle.PickleError:
    547             _log.warning("Could not pickle model, sampling singlethreaded.")

~/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/pymc3/sampling.py in _mp_sample(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, callback, discard_tuned_samples, mp_ctx, pickle_backend, **kwargs)
   1510     except KeyboardInterrupt:
   1511         if discard_tuned_samples:
-> 1512             traces, length = _choose_chains(traces, tune)
   1513         else:
   1514             traces, length = _choose_chains(traces, 0)

~/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/pymc3/sampling.py in _choose_chains(traces, tune)
   1528     lengths = [max(0, len(trace) - tune) for trace in traces]
   1529     if not sum(lengths):
-> 1530         raise ValueError("Not enough samples to build a trace.")
   1531 
   1532     idxs = np.argsort(lengths)[::-1]

ValueError: Not enough samples to build a trace.

Analysis stopped here because sampling did not converge. As the plot shows, some data points are very far away from the others, which would require the analysis to be based on more heavy-tailed distributions.

HAsfc9

with pm.Model() as model:
    HAsfc9Model = TwoFactorModel('HAsfc9',x1,x2,dataZ["HAsfc9_z"].values)

Verify model settings

HAsfc9Model.printParams(x1,x2,dataZ["HAsfc9_z"].values)

The number of levels of the x variables are (2, 11)
The standard deviations used for the beta priors are  (1.2761156280868398, 3.214359712024942)
The standard deviations used for the M12 priors are 0.17292035344730636

pm.model_to_graphviz(HAsfc9Model)

Check prior choice

with HAsfc9Model as model:
    prior_pred_HAsfc9 = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)

plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_HAsfc9,dataZ["HAsfc9_z"].values,'HAsfc9')

Prior choice is as intended: Broad over the data range.

Sampling

with HAsfc9Model as model:
    trace_HAsfc9 = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.99,random_seed=random_seed)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [HAsfc9_M12_dist, HAsfc9_mu_M12, HAsfc9_b2_dist, HAsfc9_b2_beta, HAsfc9_b1_dist, HAsfc9_b0_dist, HAsfc9_sigmaY, HAsfc9_nuY, HAsfc9_sigma12, HAsfc9_mu_b2, HAsfc9_mu_b1, HAsfc9_mu_b0, HAsfc9_sigma2, HAsfc9_sigma1, HAsfc9_sigma0]

100.00% [8000/8000 14:02<00:00 Sampling 4 chains, 0 divergences]

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 842 seconds.

with HAsfc9Model as model:
    if writeOut:
        with open(outPathData + 'model_{}.pkl'.format('HAsfc9'), 'wb') as buff:
            pickle.dump({'model': HAsfc9Model, 'trace': trace_HAsfc9}, buff)

Check sampling

with HAsfc9Model as model:
    dataTrace_HAsfc9 = az.from_pymc3(trace=trace_HAsfc9)

pm.summary(dataTrace_HAsfc9,hdi_prob=0.95).round(2)

	mean	sd	hdi_2.5%	hdi_97.5%	mcse_mean	mcse_sd	ess_mean	ess_sd	ess_bulk	ess_tail	r_hat
HAsfc9_mu_b0	-0.08	0.81	-1.71	1.52	0.01	0.01	4112.0	2349.0	4119.0	3163.0	1.0
HAsfc9_mu_b1[0]	-0.03	0.74	-1.43	1.44	0.01	0.01	3669.0	2074.0	3688.0	2541.0	1.0
HAsfc9_mu_b1[1]	-0.05	0.74	-1.44	1.43	0.01	0.01	3625.0	2253.0	3627.0	2929.0	1.0
HAsfc9_mu_b2[0]	-0.12	0.66	-1.46	1.07	0.01	0.01	3161.0	2588.0	3171.0	3161.0	1.0
HAsfc9_mu_b2[1]	-0.07	0.67	-1.42	1.22	0.01	0.01	3340.0	2428.0	3338.0	2786.0	1.0
...	...	...	...	...	...	...	...	...	...	...	...
HAsfc9_mu[225]	-0.21	0.05	-0.31	-0.11	0.00	0.00	3915.0	3915.0	3920.0	3917.0	1.0
HAsfc9_mu[226]	-0.18	0.07	-0.31	-0.05	0.00	0.00	4165.0	4160.0	4175.0	3892.0	1.0
HAsfc9_mu[227]	-0.21	0.05	-0.31	-0.11	0.00	0.00	3915.0	3915.0	3920.0	3917.0	1.0
HAsfc9_mu[228]	-0.18	0.07	-0.31	-0.05	0.00	0.00	4165.0	4160.0	4175.0	3892.0	1.0
HAsfc9_mu[229]	-0.21	0.05	-0.31	-0.11	0.00	0.00	3915.0	3915.0	3920.0	3917.0	1.0

334 rows × 11 columns

plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_HAsfc9,dataTrace_HAsfc9,'HAsfc9')

/home/bob/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/arviz/plots/backends/matplotlib/pairplot.py:216: UserWarning: rcParams['plot.max_subplots'] (40) is smaller than the number of resulting pair plots with these variables, generating only a 8x8 grid
  UserWarning,

with HAsfc9Model as model:
    plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_HAsfc9,'HAsfc9')

with HAsfc9Model as model:
    plotting_lib.pm.energyplot(trace_HAsfc9)

Posterior predictive distribution

with HAsfc9Model as model:
    posterior_pred_HAsfc9 = pm.sample_posterior_predictive(trace_HAsfc9,samples=numPredSamples,random_seed=random_seed)

/home/bob/Documents/Projekt_Neuwied/SSFA/ssfa-env/lib/python3.7/site-packages/pymc3/sampling.py:1708: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample
  "samples parameter is smaller than nchains times ndraws, some draws "

100.00% [2000/2000 00:01<00:00]

plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_HAsfc9,posterior_pred_HAsfc9,dataZ["HAsfc9_z"].values,'HAsfc9')

plotting_lib.plotLevels(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_HAsfc9,'HAsfc9',x1,x2)

plotting_lib.plotLevelsStd(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_HAsfc9,'HAsfc9',x1,x2)

Compare prior and posterior for model parameters

with HAsfc9Model as model:
    pm_data_HAsfc9 = az.from_pymc3(trace=trace_HAsfc9,prior=prior_pred_HAsfc9,posterior_predictive=posterior_pred_HAsfc9)

arviz.data.io_pymc3 - WARNING - posterior predictive variable HAsfc9_y's shape not compatible with number of chains and draws. This can mean that some draws or even whole chains are not represented.

plotting_lib.plotPriorPosteriorB(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,pm_data_HAsfc9,'HAsfc9')

Posterior and contrasts

plotting_lib.plotPosterior(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,pm_data_HAsfc9,'HAsfc9')

df_hdi_HAsfc9 = plotting_lib.plotTreatmentPosterior(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_HAsfc9,'HAsfc9',x1,x2)

df_hdi_HAsfc9

	Treatment_i	Treatment_j	hdi_ConfoMap_2.5%	hdi_ConfoMap_97.5%	isSignificant_on_ConfoMap	hdi_Toothfrax_2.5%	hdi_Toothfrax_97.5%	isSignificant_on_Toothfrax
0	Dry grass	Dry bamboo	0.033511	0.396677	True	0.027804	0.321135	True
1	Dry lucerne	Dry bamboo	-0.057178	0.067876	False	-0.059019	0.061239	False
2	Dry lucerne	Dry grass	-0.380004	-0.016990	True	-0.317744	-0.026117	True
3	BrushNoDirt	BrushDirt	-0.091347	0.223721	False	-0.115916	0.280188	False
4	Control	BrushDirt	-0.010468	0.410746	False	-0.038884	0.471882	False
5	Control	BrushNoDirt	-0.105794	0.376993	False	-0.135469	0.472578	False
6	RubDirt	BrushDirt	-0.094686	0.182720	False	-0.125451	0.167933	False
7	RubDirt	BrushNoDirt	-0.190732	0.158962	False	-0.271961	0.170081	False
8	RubDirt	Control	-0.371713	0.072724	False	-0.480866	0.068933	False
9	Clover+dust	Clover	-0.128477	0.165133	False	-0.085937	0.180708	False
10	Grass	Clover	-0.304313	0.064998	False	-0.139161	0.119079	False
11	Grass	Clover+dust	-0.276306	0.045259	False	-0.172873	0.076407	False
12	Grass+dust	Clover	-0.329568	0.009422	False	-0.225249	0.033805	False
13	Grass+dust	Clover+dust	-0.311129	-0.032240	True	-0.253405	0.000102	False
14	Grass+dust	Grass	-0.234865	0.131714	False	-0.210559	0.036078	False

if writeOut:
    df_hdi_HAsfc9.to_csv(outPathData+ 'hdi_{}.csv'.format('HAsfc9'))

plotting_lib.plotTreatmentPosteriorDiff(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_HAsfc9,'HAsfc9',x1,x2)

HAsfc81

with pm.Model() as model:
    HAsfc81Model = TwoFactorModel('HAsfc81',x1,x2,dataZ["HAsfc81_z"].values)

Verify model settings

HAsfc81Model.printParams(x1,x2,dataZ["HAsfc81_z"].values)

pm.model_to_graphviz(HAsfc81Model)

Check prior choice

with HAsfc81Model as model:
    prior_pred_HAsfc81 = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)

plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_HAsfc81,dataZ["HAsfc81_z"].values,'HAsfc81')

Prior choice is as intended: Broad over the data range.

Sampling

with HAsfc81Model as model:
    trace_HAsfc81 = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.99,random_seed=random_seed)

with HAsfc81Model as model:
    if writeOut:
        with open(outPathData + 'model_{}.pkl'.format('HAsfc81'), 'wb') as buff:
            pickle.dump({'model': HAsfc81Model, 'trace': trace_HAsfc81}, buff)

Check sampling

with HAsfc81Model as model:
    dataTrace_HAsfc81 = az.from_pymc3(trace=trace_HAsfc81)

pm.summary(dataTrace_HAsfc81,hdi_prob=0.95).round(2)

plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_HAsfc81,dataTrace_HAsfc81,'HAsfc81')

with HAsfc81Model as model:
    plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_HAsfc81,'HAsfc81')

with HAsfc81Model as model:
    plotting_lib.pm.energyplot(trace_HAsfc81)

Posterior predictive distribution

with HAsfc81Model as model:
    posterior_pred_HAsfc81 = pm.sample_posterior_predictive(trace_HAsfc81,samples=numPredSamples,random_seed=random_seed)

plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_HAsfc81,posterior_pred_HAsfc81,dataZ["HAsfc81_z"].values,'HAsfc81')

plotting_lib.plotLevels(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_HAsfc81,'HAsfc81',x1,x2)

Compare prior and posterior for model parameters

with HAsfc81Model as model:
    pm_data_HAsfc81 = az.from_pymc3(trace=trace_HAsfc81,prior=prior_pred_HAsfc81,posterior_predictive=posterior_pred_HAsfc81)

plotting_lib.plotPriorPosteriorB(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,pm_data_HAsfc81,'HAsfc81')

plotting_lib.plotLevelsStd(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_HAsfc81,'HAsfc81',x1,x2)

Posterior and contrasts

plotting_lib.plotPosterior(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,pm_data_HAsfc81,'HAsfc81')

df_hdi_HAsfc81 = plotting_lib.plotTreatmentPosterior(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,dictTreatment,dictSoftware,trace_HAsfc81,'HAsfc81',x1,x2)

df_hdi_HAsfc81

if writeOut:
    df_hdi_HAsfc81.to_csv(outPathData+ 'hdi_{}.csv'.format('HAsfc81'))

Bimodal distribution in contrast plots of HAsfc81

For e.g. the pair Clover+Dust vs. Clover shows an unexpected bimodal distribution of the contrast. We now examine the traces carefully to exclude errors in the sampling:
Get the traces on ConfoMap of the interactions

m12_clover = trace_HAsfc81['HAsfc81_M12'][:,0,2]
m12_cloverDust = trace_HAsfc81['HAsfc81_M12'][:,0,3]
diff_m12 = m12_cloverDust-m12_clover

Get the traces of the treatments:

b2_clover = trace_HAsfc81['HAsfc81_b2'][:,2]
b2_cloverDust = trace_HAsfc81['HAsfc81_b2'][:,3]
diff_b2 = b2_cloverDust-b2_clover

Look at all the pairs

sns.pairplot(data=pd.DataFrame.from_dict(
    {'m12_clover':m12_clover,'m12_cloverDust':m12_cloverDust,'diff_m12':diff_m12,
     'b2_clover':b2_clover,'b2_cloverDust':b2_cloverDust,'diff_b2':diff_b2
    }
));

Let's have that of the differences again

sns.jointplot(x=diff_b2,y=diff_m12);

We see two sets that are distributed along parallel lines in the scatter plot. This means that the model estimates two subsets of possible differences.
However, when looking at the raw data at 'R_analysis/plots/SSFA_Sheeps_plot.pdf' one can see that the distribution of values for HAsfc81 on Clover (Sheep) measured by ConfoMap appear to have a bimodal distribution.
Thus, in combination with the chosen uninformative priors, the model correctly describes the distribution as bimodal.
In summary, we see no issue with the modeling and sampling.

Summary

Set the surface parameters for every treatment dataframe:

df_hdi_Asfc["SurfaceParameter"] = "Asfc"
df_hdi_HAsfc9["SurfaceParameter"] = "HAsfc9"
#df_hdi_HAsfc81["SurfaceParameter"] = "HAsfc81"
df_hdi_R["SurfaceParameter"] = "R²"
df_hdi_epLsar["SurfaceParameter"] = "epLsar"

df_hdi_total = pd.concat([df_hdi_epLsar,df_hdi_R,df_hdi_Asfc,df_hdi_HAsfc9],ignore_index=True)

Show the treatment pairs and surface parameters where the softwares differ

df_summary = df_hdi_total[df_hdi_total.isSignificant_on_ConfoMap != df_hdi_total.isSignificant_on_Toothfrax][["Treatment_i","Treatment_j","SurfaceParameter","isSignificant_on_ConfoMap","isSignificant_on_Toothfrax","hdi_ConfoMap_2.5%","hdi_ConfoMap_97.5%","hdi_Toothfrax_2.5%","hdi_Toothfrax_97.5%"]]
df_summary

	Treatment_i	Treatment_j	SurfaceParameter	isSignificant_on_ConfoMap	isSignificant_on_Toothfrax	hdi_ConfoMap_2.5%	hdi_ConfoMap_97.5%	hdi_Toothfrax_2.5%	hdi_Toothfrax_97.5%
10	Grass	Clover	epLsar	False	True	-0.000091	0.003047	0.000486	0.003751
11	Grass	Clover+dust	epLsar	False	True	-0.000281	0.002971	0.000105	0.003410
15	Dry grass	Dry bamboo	R²	True	False	0.001541	0.002830	-0.000339	0.000455
16	Dry lucerne	Dry bamboo	R²	True	False	0.000551	0.001855	-0.000373	0.000178
17	Dry lucerne	Dry grass	R²	True	False	-0.001538	-0.000439	-0.000609	0.000202
21	RubDirt	BrushDirt	R²	True	False	-0.001555	-0.000071	-0.001123	0.000027
26	Grass	Clover+dust	R²	False	True	-0.000705	0.000535	-0.000831	-0.000136
38	RubDirt	Control	Asfc	True	False	-44.050859	-2.428273	-37.081175	0.973047
58	Grass+dust	Clover+dust	HAsfc9	True	False	-0.311129	-0.032240	-0.253405	0.000102

if writeOut:
    df_summary.to_csv(outPathData+ 'summary_filter_strong.csv')

Write out

!jupyter nbconvert --to html Statistical_Model_TwoFactor_filter_strong.ipynb

!jupyter nbconvert --to markdown Statistical_Model_TwoFactor_filter_strong.ipynb